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ABSTRACT 

Using a central difference scheme, it is necessary to add an artificial viscosity in order 
to reach a steady state. This viscosity usually consists of a linear fourth difference to 
eliminate odd-even oscillations and a nonlinear second difference to suppress oscillations 
in the neighborhood of steep gradients. There are free constants in these differences. As 
one increases the artificial viscosity, the high modes are dissipated more and the scheme 
converges more rapidly. However, this higher level of viscosity smooths the shocks and 
eliminates other features of the flow. Thus, there is a conflict between the requirements of 
accuracy and efficiency. Examples are presented for a variety of three-dimensional inviscid 
solutions over isolated wings. 
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I. Introduction 


We analyze a finite volume cell-centered scheme to solve the three dimensional Eu- 
ler equations. For a uniform Cartesian mesh the scheme reduces to a standard central 
difference scheme. Hence, one needs to add an artificial viscosity to prevent even-odd 
oscillations and also to suppress oscillations in the neighborhood of steep gradients. In 
order to accelerate the convergence to a steady state several acceleration techniques are 
used. These include, local time steps, residual smoothing and a multigrid strategy. We 
also describe other changes to the original code that either increase the accuracy of the 
steady solution or else improve the convergence rate of the iteration process. 

In order for the multigrid scheme to work it is essential that the errors be smoothed by 
the relaxation technique. Since, a central difference scheme does not include any dissipation 
one needs to add an artificial viscosity to damp the high modes. Hence, the artificial 
viscosity is needed both to give the correct shock structure in the steady state and also to 
eliminate high modes so that one can pass to a coarser mesh. As one increases the level 
of the artificial viscosity (up to some maximum) the high modes are more dissipated and 
the scheme converges more rapidly. However, this higher level of viscosity smooths the 
shocks and eliminates other features of the flow. Hence, there arises a conflict between the 
requirements of accuracy and the need to reach a steady state rapidly. 


II. Finite Volume Formulation 

The Euler equations for an inviscid compressible flow can be written in divergence form 


as 


where 


and for an ideal gas 


dQ df dg dh n 
dt dx dy dz 

(1) 

Q = (p,pu,pv,pw,E)* 

(2a) 

f = (pu, pu 2 + p, puv , puw, ( E + p)u)‘ 

(26) 

g = ( pv , puv, pv 2 + p, pvw, (E + p)v )‘ 

(2c) 

h = ( pw , puw, pvw, pw 2 + p, (E + p)w )* 

(2d) 

p=(' 1 -l)[E-p(u 2 + v 2 + w 2 )/2). 

(2e) 

the form 


^ + div(F) = 0. 

(16) 
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We integrate (1) over a three dimensional cell and consider Qij t k as an approximation 
to the average of Q over the cell. Hence, 

dQijjk , fJJdivFdV 
dt fffdV 

or using the divergence theorem, 

§~ t {VQhk + I f F-ndS = 0. (3) 

Hence, the time change of the average Q is governed by the fluxes entering and leaving the 
cell [4,5]. 

One can arrive at a similar scheme by introducing new coordinates 


£ = t{x,y,z) 

n = v(x,y,z) (4a) 

£ = c(*>y,2) 


such that £, rj, f = constant represent coordinate surfaces. Using this mapping technique 
together with finite differences leads to a formula similar to (3). However, now the volume 
is replaced by a Jacobian, 


d(x,y,z) 


(46) 


For an infinitesimally small cell the volume is equal to the Jacobian. Also, in two space 
dimensions the area of a quadrilateral is exactly equal to a central difference formula for 
the Jacobian i.e., if 


-j.J'+t ~ (&+1.J+1 - &,,+ 1 + &+i,i ~ ~ Vi+l.j + ~ / 4 

(5) 

“(&+i,/+l - Zi+I,i + Zi,j+1 - $ij) • (r}i+i,j + i - rji,)+ 1 + Vi+i,j - V i,)/4. 

Jij gives the exact area of the quadrilateral with corners (Cf.y, »?,-,/), (6+i,y»*7f+i,i)> 
(ft.j+ii Vij+i)) an d (C«+i,y+i> *7*+i,y+i)- However, in three dimensions there are differences 
between the finite volume formulation (3) and the finite difference formulation based on 
(4). In three dimensions one cannot find the volume of a three dimensional quadrilateral. 
If we assume that each face lies in a plane then we can divide the three dimensional cell into 
six pyramids and so calculate the volume. If the faces do not lie in a plane then this gives 
an approximation to the volume. This approximation is no longer the same as that given 
by central differences of the Jacobian. Similar differences occur in the flux terms where 
normals to surfaces are required. The formula that comes from a finite volume approach is 
no longer exactly the same as that of a finite difference plus mapping approach. However, 
both formulas agree to second order accuracy. In a Cartesian mesh the finite volume and 
cell-centered finite difference approaches are the same. 
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III. Artificial Viscosity 


Both the finite volume and the finite difference approaches lead to a pure central 
difference method for Cartesian grids. Though this scheme is stable for constant coefficient 
hyperbolic equations it is subject to instabilities that will prevent the convergence to 
a steady state. To force this convergence a fourth difference viscosity is added to the 
scheme. The fourth difference causes oscillations in the neighborhood of shocks. Hence, a 
nonlinear second difference is added to control oscillations near the shocks and the fourth 
difference is turned off. The total artificial viscosity, V, is the sum of such second and 
fourth differences in each coordinate direction. 


V M 







_l_y* , _ v* i 

' V *,},k+ h v k 


( 6 ) 


Hence it is sufficient to describe these terms in the £ direction. Since we only take differ- 
ences at neighboring points the artificial viscosity is always in conservation form. 

The first difference is defined as 


Di+lj'k — Qi+l,j,k Qi,),k 


(7a) 


and the second £ difference is defined as 

Ei,},* = “ A- $,/,*• ( 7 &) 

We then form the second and fourth differences. In particular the fourth difference is 
formed as a second difference of a second difference with positive weights [3,8]. Hence, 




(<) 


cW 


( 8 ) 


Let, 




Pj+l, j,k 2Pj,j,k Pj—l,j,k 


| J’l+l.i.t + + Pl-l,j,k | <9a) 

Then is used to detect the location of shocks. When is large then the fourth 
difference is reduced. Other ways of normalizing the second difference of the pressure are 
also possible. Let, 


a i+y,k = k{2) I'i+ij*, ^+2,i.fc) • 


(96) 


We also multiply a by a function of the Mach number to reduce a near the surface. Finally 
let A be a measure of the fluxes (this will be discussed in more detail). Then 


6. 


( 2 ) 


i+h,j,k 


<7... l 


(9c) 
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4y,fc = max (o> ^ (4) - <n,i,k) • ( 9d ) 

Let, A = |£,J3 = |^,C = 1^-, where F,G,H are the fluxes in the coordinate system 
(f , rj, f). The original code chose A as 

X ( = X n = X s — p(A) + p{B) + p(C) (10a) 

where p is the the spectral radius of the matrix. For problems with a highly stretched 
mesh it was found [1,3,8,11,12] that for increased accuracy one should choose 

A* = p(A), A’ = p(B), A< = p(C). (106) 

K ( J ) ) K ( 4 ) are constants that determine the level of the second and fourth differences. These 
constants are given as input to the code. In the result section these will be varied to see 
their effect both on the accuracy of the solution and on their convergence rate. 

In all our calculations we use a C grid in each plane. It was found important to calculate 
the artificial viscosity across the wake and not to treat the wake as a solid surface. In the 
spanwise direction we use either an H or an O grid. The H grid leads to a simpler topology 
but does not give sufficient resolution near the wing tip. With the O mesh care must be 
taken to reflect the proper points across the wake and near the wing tip. This is especially 
important for the fourth difference in the £ direction which requires several points on the 
other side of the wake. 


IV. Results 


We first consider flow past an isolated ONERA M6 wing with M„ 0 = 0.84 and a = 3.06°. 
We use a coarse 96 x 16 x 16 C-H mesh. In Figure la, lb we show the C p plots over the 
upper surface. For this case 

*( 2 ) = 2 . 12 . . , 

= 2./64. (llo) 

In Figure lc we plot the convergence history for this case. A five stage Runge-Kutta 
scheme was used with two evaluations of the artificial viscosity. This scheme was used once 
per iteration on the finest grid, twice on the next finest grid and three times on all coarser 
grids. On the way up from coarser to finer grids no smoothing was used and the changes 
were just interpolated to the next finer grid. With this scheme we were able to reduce the 
mean density residual by 9 orders of magnitude after 100 iterations on the finest mesh. 
Using a FMG method the initial condition was obtained by using 10 iterations each on two 
coarser grids. These convergence rates are close to those obtained by Jameson [6] using a 
nodal scheme. Hence, we see that there are no major differences in convergence rates and 
robustness between the cell-centered and the nodal versions of the multi-stage methods. 
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In Figure 2 we plot the same case but with 

/cW = 0.5/2. 

= 0.5/64. 

We see that now the shocks are much sharper and that other features of the flow are more 
pronounced. In addition, the total lift was lowered from .299 to .287 and the drag was 
changed from .0174 to 0.147. With this reduced k 1 2 \ k/ 4 ) the lift and drag are close to 
that given in the finer mesh of Figure 3. However, now the convergence rate is reduced to 
a 5 order reduction within the 100 iterations. 

In Figure 3 we plot the same physical case but using a finer 192 x 32 x 32 C-0 mesh. 
Looking carefully at individual stations one can again see that the lower viscosity level 
leads to sharper profiles but at the expense of a slightly reduced convergence rate. With 
this finer mesh there is much less of a sensitivity to the constants in the artificial viscosity 
compared with the coarse mesh (see also [12]). 

In Figure 4 we plot the convergence rate for a supersonic flow about a delta wing. 
For this case a 128 x 24 x 16 C-H mesh was used supplied by Moitra [7]. This further 
demonstrates the robustness of the present code over different configurations and different 
flight conditions. 
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Figure lb. Upper surface pressure plot. 
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Figure lc. Convergence plot of average p residual. 
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Figure 3a. Same as Figure 1 with 192 x 32 x 32 C-0 mesh. /J 2 ) = i «( 4 ) = i 

2 7 64 
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Figure 3b. Upper surface pressure. 
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Figure 3d. Convergence history for Figure 3c. Note: Scale is different than Figure lc and 
2c since fewer iterations were done. 
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Figure 4a. Delta wing = 2.5, a = 1.2°. 



WORK 


CASE 1 NMESH-3 NCYC-IOO 
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Figure 4b. Convergence history. 
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